##This file assess the robustness by varying the bandwidth within which the RD is estimated 

rm(list=ls(all=TRUE))
graphics.off()
library(rdd)
library(xtable)

setwd() # set working directory
load("data/data_alt_boot.rdata")
#drop all from closed and split lists 
data_alt <- data_alt[which(data_alt$marg.los!=0 & data_alt$marg.win != 1 & 
                             data_alt$openlist == 1 & abs(data_alt$margsim) < 1
                           & data_alt$splitlist == 0 & 
                             data_alt$simscore > 0 & data_alt$simscore < 1),]

#if cases are extremely close they might have wrong signed simscore. Recode these to arbitrarily small simscore
data_alt$margsim[data_alt$margsim <= 0 & data_alt$elected == 1] <-  0.001
data_alt$margsim[data_alt$margsim >= 0 & data_alt$elected == 0] <- -0.001

data_alt <- na.omit(data.frame(data_alt[c("margsim", "electedt1", "rerun")]))

#find Imbens-Kalyanaraman bandwidth
IK.bw  <- IKbandwidth(data_alt$margsim, data_alt$electedt1)
IK.bw2 <- IKbandwidth(data_alt$margsim, data_alt$rerun)

##Evaluate robustness by varying bandwidth 

rollests1    <- rep(NA, ceiling(100*IK.bw*2))[-(1:4)]
rollses1     <- rep(NA, ceiling(100*IK.bw*2))[-(1:4)]
rollbounds1  <- cbind(rollests1,rollests1)

Ns           <- rep(NA, ceiling(100*IK.bw2*2))[-(1:4)]

for (i in 4+1:length(rollests1)){
  rd                 <- RDestimate(electedt1 ~ margsim, data = data_alt, bw = i/100)
  rollests1[i-4]     <- rd$est[1]
  rollbounds1[i-4,]  <- c(as.numeric(rollests1[i-4]) + rd$se[1]*1.96,
                          as.numeric(rollests1[i-4]) - rd$se[1]*1.96)
  rollses1[i-4]      <- rd$se[1]
  Ns[i-4]            <- rd$obs[1]
}


quartz(w=8, h =4)
par(mfrow=c(1,2), 
    mar=c(4,4,2,1),
    cex = 1.1)
plot(1:length(rollests1)/100+0.04, rollests1, ylim=c(-.5,.5),pch = 19,
     cex = 0.1,
     xlab="Bandwidth",ylab="Incumbency advantage", 
     cex.main = 0.9,
     main = "Effect on p(rerunning and elected)")
lines(1:length(rollests1)/100+0.04, rollests1)
lines(1:length(rollests1)/100+0.04, rollbounds1[,1])
lines(1:length(rollests1)/100+0.04, rollbounds1[,2])
abline(h=0, col="grey50")
abline(v=IK.bw/2, col="grey50", lty = 2)
abline(v=IK.bw, col="grey50", lty = 3)
abline(v=IK.bw*2, col="grey50", lty = 4)

legend("bottomright", lty = c(2:4), c("half-IK-BW", "IK-BW", "2*IK-BW"), 
       col=c("grey50","grey50","grey50"), box.lwd = 0, box.col = "white", bg = "white")
box()

#Repeat for rerunning

rollests2    <- rep(NA, ceiling(100*IK.bw2*2))[-(1:4)]
rollses2     <- rep(NA, ceiling(100*IK.bw2*2))[-(1:4)]
rollbounds2  <- cbind(rollests2,rollests2)

for (i in 4+1:length(rollests2)){
  rd                 <- RDestimate(rerun ~ margsim, data = data_alt, bw = i/100)
  rollests2[i-4]     <- rd$est[1]
  rollbounds2[i-4,]  <- c(as.numeric(rollests2[i-4]) + rd$se[1]*1.96,
                          as.numeric(rollests2[i-4]) - rd$se[1]*1.96)
  rollses2[i-4]      <- rd$se[1]
}

plot(1:length(rollests2)/100+0.04, rollests2, ylim=c(-.5,.5),pch = 19,
     cex = 0.1,
     xlab="Bandwidth",ylab="Incumbency advantage",
     cex.main = 0.9,
     main = "Effect on p(rerunning)")
lines(1:length(rollests2)/100+0.04, rollbounds2[,1])
lines(1:length(rollests2)/100+0.04, rollbounds2[,2])
abline(h=0, col="grey50")
abline(v=IK.bw2/2, col="grey50", lty = 2)
abline(v=IK.bw2, col="grey50", lty = 3)
abline(v=IK.bw2*2, col="grey50", lty = 4)

legend("bottomright", lty = c(2:4), c("half-IK-BW", "IK-BW", "2*IK-BW"), 
       col=c("grey50","grey50","grey50"), box.lwd = 0, box.col = "white", bg = "white")
box()

